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Water  management  has  long  been  recognized  as  a  critical  issue  in  the  operation  of  a  Polymer  Electrolyte 
Membrane  Fuel  Cell  (PEMFC).  If  the  membrane  is  allowed  to  dehydrate  then  ionic  conductivity  will  drop 
and  result  in  significant  power  losses.  At  the  opposite  extreme,  if  not  enough  water  is  removed  from  the 
membrane  then  liquid  water  will  accumulate  in  the  Gas  Diffusion  Layer  (GDL)  and  block  the  transport 
of  reactants  to  the  reaction  sites.  In  this  work,  the  dynamics  of  membrane  hydration  are  analyzed  while 
under  the  influence  of  a  closed-loop  control  system.  From  a  controller  assessment  perspective,  the  mem¬ 
brane  model  is  unique  in  that  through  the  plane  spacially  dependent  water  accumulation  is  captured. 
By  combining  with  an  electrochemical  model  and  simple  material  and  energy  balances  over  the  solid 
and  fluid  materials,  the  dynamics  of  the  membrane  are  shown  to  be  strongly  influenced  by  the  thermal 
responses  of  the  solid  as  well  as  humidity  levels  in  the  gas  streams,  all  of  which  are  a  function  of  the 
controller  utilized.  We  conclude  by  illustrating  the  highly  sensitive  nature  of  the  system  with  respect  to 
water  diffusivity  within  the  membrane. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  objective  of  a  Polymer  Electrolyte  Membrane  Fuel  Cell 
(PEMFC)  is  to  convert  hydrogen  into  electric  power.  Central  to 
PEMFC  efficiency  is  the  ionic  conductivity  of  the  membrane,  which 
is  strongly  influenced  by  membrane  hydration  levels.  Specifically, 
greater  hydration  will  result  in  greater  conductivity  and  thus  a 
more  efficient  cell.  However,  if  hydration  levels  exceed  the  capacity 
of  the  membrane,  then  a  layer  of  liquid  water  will  begin  to  form  and 
block  the  transfer  of  reactants  to  the  reaction  sites.  Thus,  the  oper¬ 
ational  objective  with  regard  to  membrane  hydration  is  to  operate 
at  a  level  just  below  the  flooding  limit.  Unfortunately,  the  antici¬ 
pated  applications  of  a  PEMFC  (most  notably  automotive)  suggest 
that  frequent  changes  in  power  demand  will  be  the  norm.  Since 
the  reaction  product  is  water,  rate  changes  in  power  demand  will 
result  in  changes  in  water  production  within  the  membrane. 

From  a  steady-state  perspective,  the  water  management  chal¬ 
lenge  is  to  ensure  that  water  removal  rates  are  equal  to  production. 
While  water  production  is  proportional  to  current  density,  the  flux 
of  water  from  the  membrane  is  a  complicated  function  of  operating 
temperature  and  hydration  level  (within  the  membrane  as  well  as 
within  the  reactant  gases).  Within  the  dynamic  framework,  there  is 
the  additional  question  of  hydration  level  excursions  during  tran¬ 
sient  periods.  In  some  cases,  the  system  may  not  recover  from  an 
excursion  into  flooding  or  dehydration,  even  if  a  steady-state  based 


*  Corresponding  author.  Tel.:  +1  312  567  3537;  fax:  +1  312  567  8874. 
E-mail  address:  chmielewski@iit.edu  (D.J.  Chmielewski). 

0378-7753 /$  -  see  front  matter  ©  2011  Elsevier  B.V.  All  rights  reserved, 
doi:  1 0.1 01 6/j.jpowsour.201 1 .01 .096 


analysis  suggests  otherwise.  Within  these  transient  scenarios  key 
factors  include  the  hydration  capacity  of  the  membrane  and  the 
flowing  gases  as  well  as  thermal  response  time. 

Due  to  the  diversity  of  phenomena  occurring  within  a  rather 
complex  structure,  a  wide  variety  of  PEMFC  models  can  be  devel¬ 
oped,  each  focusing  on  a  different  dimension  and  time/length  scale. 
Overviews  of  PEMFC  modeling  can  be  found  in  [1-4].  In  the  fol¬ 
lowing  we  will  focus  on  models  of  the  membrane.  Initial  efforts  to 
model  the  membrane  in  a  PEMFC  include  Springer  et  al.  [5]  and 
Bernardi  and  Verbrugge  [6].  Concerning  flux  of  water  through  the 
membrane  both  include  the  electro-osmotic  drag  term,  which  is 
proportional  to  current  density  and  causes  water  to  move  toward 
the  cathode.  In  the  Springer  et  al.  [5]  model,  the  flux  also  includes 
a  back  diffusion  term,  which  serves  to  counter  act  the  impact 
of  electro-osmotic  drag.  In  the  Bernardi  and  Verbrugge  model,  a 
hydraulic  pressure  term  is  used  to  counteract  electro-osmotic  drag. 
Efforts  that  directly  employ  the  Bernardi  and  Verbrugge  membrane 
model  include  Eikerling  et  al.  [7]  andBaschukand  Li  [8].  The  work  of 
Springer  et  al.  [5]  also  proposed  hydration  dependent  relationships 
for  the  electro-osmotic  drag  coefficient,  the  diffusion  coefficient 
and  the  ionic  conductivity  of  the  membrane. 

Nguyen  and  White  [9]  extended  the  one-dimensional  (through 
the  membrane)  model  of  [5]  to  include  an  along  the  channel  dimen¬ 
sion.  However,  to  reduce  model  complexity  the  differential  aspects 
of  the  through  the  membrane  direction  were  approximated  by 
algebraic  relations.  Specifically,  diffusion  flux  became  a  function 
of  the  hydration  difference  between  cathode  and  anode,  and  ionic 
conductivity  became  a  function  of  the  hydration  average  between 
anode  and  cathode.  The  effort  by  Yi  and  Nguyen  [10]  extended  the 
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Nomenclature 

Super-  and  sub-scripts 

a ,  cj ,  e ,  s,  m  anode,  cathode,  jacket,  ambient,  solids,  mem¬ 
brane 

J  flux  (mol  cm-2  s-1) 

C  concentration  (mol  cm-3 ) 

T  temperature  (I<) 

F  volumetric  flow  (cm3  s-1 ) 

r  reaction  rate 

j  current  density  (A  cm-2 ) 

E  voltage  (V) 

7 Z  resistance  of  membrane  ( ficm2 ) 

a  conductivity  in  membrane  (S/m2) 

Pe  power  density  ( W  cm-2 ) 

aw  activity  of  water  in  the  membrane 

P  Partial  pressure  (atm) 


model  of  [9]  by  including  the  hydraulic  pressure  term  (of  Bernardi 
and  Verbrugge  [6])  in  the  water  flux  expression.  Again  the  differ¬ 
ential  aspect  of  the  through  the  membrane  direction  was  replaced 
by  algebraic  relations.  Efforts  that  employ  the  Yi  and  Nguyen  [10] 
approach  include  Rowe  and  Li  [11],  Wu  et  al.  [12],  and  Zhou  et  al. 
[13].  In  You  and  Liu  [14]  the  membrane  model  and  cell  configura¬ 
tion  of  Yi  and  Nguyen  [10]  is  revisited,  but  the  differential  aspects 
of  the  through  the  membrane  direction  are  retained  and  solved  via 
a  Computation  Fluid  Dynamics  (CFD)  approach. 

All  of  the  models  discussed  up  to  this  point  have  been  of  the 
steady-state  variety.  The  first  class  of  dynamic  models  is  those  that 
place  dynamic  descriptions  of  temperature  and  gas  composition 
on  top  of  a  steady-state  model  of  the  membrane.  For  the  most  part, 
papers  in  this  class  assume  algebraic  relations  for  the  through  the 
plane  direction  and  view  the  cell  as  being  a  single  lump  (Pukrush- 
pan  et  al.  [  1 5  ],  Lauzze  and  Chmielewski  [16],  and  Zhang  et  al.  [  1 7  ] )  or 
as  having  an  along  the  channel  spacial  direction  (Golbert  and  Lewin 
[18,19]  and  Methekar  et  al.  [20]).  In  De  Francesco  and  Arate  [21]  a 
differential  perspective  is  used  in  the  through  the  plane  direction 
while  a  single  lump  is  employed  for  the  along  the  channel  direction. 

The  second  class  of  dynamic  models  are  those  that  try  to  capture 
the  membrane’s  ability  to  store  water  and  quantify  the  time  rate  of 
change  of  this  stored  water.  In  Shan  and  Choe  [22]  the  entire  mem¬ 
brane  is  considered  a  single  lump  and  water  flux  terms  at  the  anode 
and  cathode  sides  are  used  to  quantify  the  accumulation  of  water. 
In  Chia  et  al.  [23],  a  similar  single  lump  configuration  is  used,  and 
then  extended  to  a  lumps  in  series  configuration  to  approximate 
the  along  the  channel  spacial  direction.  In  Chen  et  al.  [24],  water 
accumulation  in  the  through  the  plane  direction  is  quantified.  In 
Wang  and  Wang  [25]  and  Um  and  Wang  [26],  CFD  methods  are 
used  to  simulate  water  accumulation  in  both  the  through  the  plane 
direction  and  the  along  the  channel  direction. 

Concerning  experimental  studies  of  membrane  hydration, 
Bellows  et  al.  [27],  have  used  neutron  imagining  methods  to 
measure  hydration  profiles  in  the  through  the  plane  direction. 
Unfortunately,  limitations  in  measurement  fidelity  have  made  it 
challenging  to  characterize  hydration  dynamics  in  the  through  the 
plane  direction  [28,29]. 

The  objective  of  this  paper  is  to  develop  an  accumulation  based 
through  the  plane  membrane  model,  similar  to  those  of  [24-26], 
but  then  analyze  this  model  within  a  closed-loop  perspective.  As 
illustrated  by  the  above  literature  review,  closed-loop  analysis  of 
an  accumulation  type  model  seems  to  be  missing  from  the  PEMFC 
modeling  as  well  as  the  PEMFC  control  literature.  The  current  effort 
aims  to  fill  this  gap. 


Cooling  Jacket 


The  paper  is  organized  as  follows.  The  PEMFC  model  is  presented 
in  Section  2,  and  it  is  solutions  procedure  is  discussed  in  Section  3. 
Analysis  of  the  model  under  a  variety  of  closed-loop  control  con¬ 
figurations  is  given  in  Section  4.  Section  5  presents  some  additional 
discussion  concerning  the  topics  of  flooding  and  diffusion  within 
the  membrane. 

2.  PEMFC  model 

The  system  scenario  is  similar  to  that  of  Lauzze  and  Chmielewski 
[16].  From  a  global  perspective  the  PEMFC  stack  is  assumed  to  be  of 
sufficient  size  that  air  cooling  is  required  (approximately  10kWe). 
However,  the  model  presented  reflects  the  volume  and  surface 
areas  of  a  single  flow  channel,  under  the  assumption  that  macro¬ 
scopic  stack  values  for  power,  current,  and  flow  rate  can  be  obtained 
by  appropriate  multiplication  of  this  modeling  unit.  Additionally, 
the  spacial  aspect  of  the  single  flow  channel  will  be  ignored  in  favor 
of  the  simplicity  of  a  Continuous  Stirred  Tank  Reactor  (CSTR)  form. 
In  spite  of  this  neglect  of  along  the  channel  spacial  dependance,  we 
have  found  this  model  to  exhibit  sufficient  richness.  In  contrast  to 
[16],  the  new  model  considers  an  open-ended  humidified  hydrogen 
feed  and  of  course  hydration  dynamics  within  the  membrane. 

The  unit  cell  of  the  model  consist  of  two  gas  chambers  separated 
by  a  membrane  electrode  assembly  (MEA),  see  Fig.  1.  On  the  anode 
side,  hydrogen  is  split  into  hydrogen  ions  and  electrons.  While  the 
ions  travel  through  the  membrane,  the  electrons  travel  through 
the  catalyst  layer  and  the  Gas  Diffusion  Layer  (GDL)  to  the  current 
collector  and  on  to  the  load.  These  electrons  then  travel  back  to 
the  cathode  where  they  combine  with  the  hydrogen  ions  and  oxy¬ 
gen  to  produce  water.  The  rate  of  reaction  is  proportional  to  the 
current  density  j/nJr=  -rH2  =  -1  /2 r0l  =  rH2o  where  r,  represents 
the  generation  of  species  i  per  unit  area  of  membrane.  While  the 
membrane  is  designed  to  be  impermeable  to  H2  and  02,  it  is  capa¬ 
ble  of  significant  water  uptake.  As  such  rH2o  cannot  be  used  for  the 
gas  phase  material  balances.  Instead,  a  pair  of  water  transfer  fluxes 
to  the  membrane  from  the  anode  and  cathode  gas  chambers,  J^0 
and 7h2o’  will  be  defined. 

While  flowrate  through  the  cathode  is  a  good  heat  removal 
mechanism,  it  is  common  to  augment  the  system  with  a  cooling 
jacket,  the  third  chamber  of  Fig.  1.  On  the  subject  of  heat  trans¬ 
fer,  the  surface  area  for  heat  transfer  in  the  cathode  and  anode 
chambers  is  more  than  just  the  membrane  surface  area,  Am.  Specifi¬ 
cally,  the  current  collector  flow  channel  walls  that  surround  the  gas 
chambers  are  also  available  for  heat  transfer.  As  such,  we  assume  all 
of  this  solid  material  including  the  membrane  to  be  a  single  lump  for 
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energy  balance  purposes.  Finally,  heat  transfer  can  also  occur  with 
the  environment,  specifically  losses  from  the  stack  edges.  Thus,  an 
effective  surface  area  is  assumed  based  on  the  expected  ratio  of 
insulated  surface  area  to  the  stack  volume,  and  then  applied  to  the 
volume  of  the  unit  cell. 

2.1.  Material  and  energy  balances 


The  heat  generation  term  Qgen  is  the  amount  of  heat  produced  by 
the  electrochemical  reaction,  given  by  Qgen  =  (AHjH2o)rH2o  -Pe. 
where  Pe  =jEceu .  It  is  additionally  noted  that  the  above  balances 
assume  positive  values  for  J^Q  and  7h2o*  ^  either  flux  is  negative 
then  the  appropriate  terms  are  replaced  by J^2qTs  or Jh2o^S  (3)’ 
(7)  and  (10). 


Using  the  above  description,  the  following  material  and  energy 
balances  are  developed.  In  the  anode  chamber: 


Vo^=F“c“20-F“c“2+rH2A„ 

dCa 

\/  h2°  _  para  pa^a  ja  a 

Va  dt  ~  to  Lh20,o  -  M  LH20  _  JH2Onm 

w  dTa  FoJo  ~  ^Ja  +  ((UA)a/CigCp4g)(Ts  -  Ta)  +  (m2Fa -J^0Ta)Am 
Va  dt  Cis 


,_fo+(% -J£2  0)An 


*1  = 


-lg 


In  the  cathode  chamber: 


dCc0 
Vc-°2 


dt  =  FoC02,o  -  F1C02  +roA 


dCf,  n 

1/  M2U  _  pcrc  Fcrc  jc  A 

vc  Atr  -roLH20,o  _MLH20  ~Jh2 0^n 


dt 


dTc  FC0TC0  -  F\TC  +  ((UA)c/CigCp>ig)(Ts  -  Tc)  +  ( r02Tc -JcH20Tc)Am 

Vc~dF  ~  a 


n= 


Fq  +  (r02  -JS2 0V»; 


At  the  solid  material  and  cooling  jacket: 

Vj  ffr  =  P'ri  -  PP  +  ^rf(Ts  -  P) 
at  (Ptp)j 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

O) 


(pC„)sVs^l  =  ( UA)a(Ta  -  P)  +  (UA)C(TC  -  P)  +  m)j(Tj  -  P) 


+  (UA)e(P  -  P)  -  (rH2TQ  +  r02P  -JS2ot“ 

— Jh20  TC)AmCpjg  +  QgenAm  (10) 


2.2.  Electrochemical  model 

The  cell  voltage  is  the  ideal  minus  losses 

Ecell  —  Ener  ~  Eact  —  —  Emt  ( 1 1 ) 

Ener  =  E°  +  (RT^/nJ=)\n{Pn2PQ^/Pii2o)  is  the  Nernst  poten¬ 
tial.  The  activation  loss  is  Eact  =  {\/a){RT^/nJ:)\n{j/j0)f 
where  j0  is  the  exchange  current  density.  The 
ohmic  loss  is  Eohm  =j7l,  where  JZ  =  fjm  dz/ar(z), 
cr(z)  =  0.005 1 93 X(z)  -  0.00326exp  ( 1 269.0(  1/303-1  /!)), 

A(z)  =  C^20(z)/Ns  and  CJ”  0(z)  is  the  hydration  level  within  the 
membrane  (defined  in  the  next  sub-section).  The  membrane  thick¬ 
ness  is  r m  {z= 0  is  anode  side  and  z  =  rm  is  the  cathode  side).  The 
mass  transfer  loss  is  Emt  =  (1/2  +  1/a)  {RP^/nP)  lnOY/Oi  -/)). 
where  ji  =  2nTkcgdlCc0^  is  the  limiting  current  density.  The  mass 

transfer  coefficient  across  the  GDL  is  kldl  =  D^/tj  where  x\  is  the 
thickness  and  Dlgdl  is  the  diffusivity  of  the  GDL  (i  =  a  or  c). 


2.3.  Membrane  hydration  model 


A  water  balance  within  the  membrane  yields: 
dC™  n  dj ™  n 

H2Q  _  _  jH2Q 

dt  dz  1  j 

where  C^0  is  the  concentration  of  water  in  the  membrane  and  JJJ  0 
is  the  flux  of  water  within  the  membrane.  Water  transport  within 
the  membrane  is  due  to  two  separate  mechanisms  -  diffusion  and 
electro-osmotic  drag: 


~  n 
J h2o  -  ~Ur 


aqy2o 

1  dz 


+  ? 


)_ 

P 


(13) 
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Fig.  3.  Typical  hydration  profiles. 


If  the  diffusion  and  drag  coefficients  Dm  and  £  are  assumed  constant, 
then  the  following  model  will  arise. 


9Cff2o 

dt 


=  D„ 


d2C! 


H20 


dz 2 


(14) 


Table  1 

Model  parameters 


Parameter 

Description,  value,  units 

T 

Faraday’s  constant,  95485  C  mol-1 

°agdi 

Gas  diffusion  in  anode  GDL,  0.1490  cm2  s_1 

°Cgdi 

Gas  diffusion  in  cathode  GDL,  0.0295  cm2  s-1 

Dm 

Water  diffusion  in  membrane,  1.5  x  10-6  cm2  s-1 

£ 

Electro-osmotic  drag  coefficient,  1 

Ns 

Number  of  sulfonic  sites,  0.00197  mol  cm-3 

Qg 

Ideal  gas  concentration,  3.45  x  10-5  molcirr3 

A 

Density  of  solids,  0.35  g  cm-3 

Cp,ig 

Ideal  gas  heat  capacity,  15  J  mol-1  K-1 

Cp,s 

Solids  heat  capacity,  0.93  J  mol-1 1<-1 

AHf;  h2o 

Heat  of  formation  of  water,  -286,  000J  mol-1 

n 

No.  of  electrons  transferred  in  reaction,  2 

a 

Charge  transfer  coefficient,  0.5 

jo 

Exchange  current  density,  0.01  mAcm-2 

E° 

Reversible  voltage,  1.2  V 

Ta 

Thickness  of  anode  GDL,  0.0350  cm 

Tc 

Thickness  of  cathode  GDL,  0.0350  cm 

rm 

Thickness  of  membrane,  0.0150  cm 

Va 

Anode  gas  volume,  1.25  cm3 

Vc 

Cathode  gas  volume,  1.25  cm3 

Vj 

Jacket  gas  volume,  7.5  cm3 

Vs 

Solid  volume,  7.5  cm3 

Am 

Membrane  area,  25  cm2 

Da,  Uc,  Uj 

Heat  transfer  coef  (gas-solid),  5.2  x  10_4J  s_1  cnrr2  K1 

Ue 

Heat  transfer  coef  (solid-ambient),  5.2  x  10_6J  s_1  cm-2  K-1 

ja 

J  H20 


+  Dn 


dC™ 


h2o 


dz 


at  z  =  0 


(15)  3.1.  Spacial  discretization 


dCm  ’t 

Dm  — ^ — b  J—  +7h2o  +  rH2o  =  0  at  z  = 


(16) 


The  flux  of  water  entering  the  membrane  from  the  gas  chambers 
(see  Fig.  2)  is  defined  as: 

Jh2o  =  J<gdl  [ch2o  -  ch2o,s] 

Jh2o  =  kgdi  [ch2o  -  ch2o,s] 

where  Cfj  0s  =  alwPvap(Ts')/RTs  and  alw  satisfies  the  gas/membrane 
equilibrium  relation  at  the  chamber  interfaces. 


The  time  dependent  nature  of  the  boundary  conditions  (15)  and 
(16)  makes  them  ill-suited  for  the  following  discretization  method. 
As  such  we  convert  (14)-(16)  into  a  homogeneous  form: 


9ch2o 


(17) 

dt 

(18) 

dCM20 

dz 

=  D„ 


3  2C\ 


dz2 


^£+/(t,Z) 


=  0  at  z  =  0  and  z  =  r+ 


(21) 

(22) 


whereat,  z)  is  selected  such  that  (21 )  and  (22)  is  equivalent  to  (14) 
and  (16).  Clearly /(z,  t)  will  need  to  be  of  the  form  ([30]): 


Ns  (o.043  +  17.81a^  -  39.85(a^)2  +  36.0(a^,)3)  =  C™0  \z=Q  (19) 

Ns  (o.043  +  17.81a^-39.85(a^)2+36.0(a^)3)  =  C™  0 1  (20) 

Fig.  3  illustrates  typical  membrane  hydration  profiles 
(Tables  1  and  2),  all  at  steady  state  and  a  solid  temper¬ 
ature  of  80  °C,  recall  that  A  =  C™20/Ns.  At  low  current 
density,  diffusion  dominates  and  results  in  a  nearly  hor¬ 
izontal  profile.  However,  at  higher  current,  the  density 
combined  effect  of  water  generation  on  the  cathode  side  along 
with  electro-osmotic  drag,  also  toward  the  cathode  side,  is 
observed. 

3.  Solution  methodology 


f(t,  z)  =  8{z)Jo{t)  +  S(z  -  rm)jzm  (23) 


Table  2 

Nominal  operating  conditions 
Parameters  Description,  value,  units 

C“2  o  Inlet  hydrogen  concentration  1.78  x  10-5  mol  cm-3 

C^Qo  Inlet  water  concentration  in  the  anode  1.66  x  10-5  mol  cm-3 

C£2  0  Inlet  oxygen  concentration  6.72  x  10-6  mol  cnrr3 

CT  n  Inlet  water  concentration  in  the  cathode  2.47  x  10-6  mol  cm-3 

H20,0 

C^2  o  Inlet  nitrogen  concentration  2.53  x  10-5  mol  cm-3 

T°  Inlet  anode  temperature  80°C 

T£  Inlet  cathode  temperature  40° C 

FJ  Initial  flowrate  of  anode  2.2  cm3  s_1 

F£  Initial  flowrate  of  cathode  5.8  cm3  s-1 

F>  Nominal  flowrate  of  the  cooling  jacket  95.9  cm3  s-1 


The  above  PEMFC  model  contains  eight  Ordinary  Differ¬ 
ential  Equations  (ODEs)  and  one  Partial  Differential  Equa¬ 
tion  (PDE).  While  finite  element  methods  could  be  used  to 
address  the  PDE  portion  of  the  model,  an  alternative  is  to 
approximate  the  PDE  by  a  set  of  ODEs.  This  approach  is 
expected  to  be  of  greater  utility  in  future  studies  aimed  at 
model  reduction  and  the  development  of  a  model  based  con¬ 
troller. 


Fig.  4.  Current  controller  configuration. 
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Fig.  5.  Simulation  of  current  controller. 


where  5(  • )  is  the  Dirac  delta  function  andj0,  ]im  are  determined  by 
integrating  (21)  fromz  =  0_  to  0+  andz  =  r~  to  which  yields: 


]o(t)  = 


Ia 

1h2o 


Ji 

T 


(24) 


irm(t)  - 


+ 


ji  +  J- 

T  nT 


(25) 


Next,  the  system  (21)  and  (22)  is  spacially  discretized  using  the 
standard  Galerkin  approach  [31  ].  Let  C^Q{t,  z )  be  approximated  as 
follows 

N 

Cff2 0(t>  z)  -  Cff2 0,N(t,  z)  =  Y)mj(t)pj(z)  (26) 

j=i 

where  pj(z)  =  Hj  cos (Wj-z)  is  a  sequence  of  basis  functions  (wj  = 

7i(j  -  l)/rm,  Hj  =  1/  jxfn  if  j  =  l  and  Hj  =  l/^/(rm/2)  otherwise). 
As  the  eigenfunction  of  (21),  these  are  known  to  be  orthonormal 
under  the  inner  product 

f‘Tm 

(Pi,  Pj )  =  /  Pi{z)Pj(z)dz  (27) 

Jo 

Now  define  a  residual  function 

f)cm  d2cm 

%  =  -^+Dm^^+/(t,z)  (28) 

and  enforce  the  conditions  (%,  pd  =  0,  i  =  1, . . .,  N.  This  results  in  the 
following  set  of  ODEs,  which  will  be  used  to  approximate  C^Q{t,  z ) 


with  the  help  of  Eq.  (26). 

^  =  -Dmwfnij  +  Pi(O)Ut)  +  p,(rm)]Tm{t)  i  =  1,...,N  (29) 


3.2.  Model  structure 


The  proposed  PEMFC  model  contains  8  +  N  ODEs 
(Eqs.  ( 1 )— (3 ),  (5)— (7),  (9),  (10),  and  (29))  along  with 

a  number  of  algebraic  relations.  In  most  cases,  these 
algebraic  relations  are  simple  functions  of  the  state 
variables: 

x  =  n2  ch2o  Ta  C£2  c^2o  r  V  r  mf  (30) 

or  the  manipulated  variables 

u  =  [F0fl  Fg  Fj  Ecell]T  (31) 

and  can  be  directly  substituted  into  the  differential  equations.  How¬ 
ever  algebraic  relations  ( 1 1 ),  ( 1 9)  and  (20)  are  such  that  an  analytic 
expression  for  the  current  density,  j,  and  the  membrane  surface 
activities,  acw  and  a^,  are  not  easily  obtained.  For  these  three  rela¬ 
tions,  a  bisection  search  algorithm  is  employed  at  each  time-step 
of  the  numeric  integration  scheme.  In  sum,  the  PEMFC  model  has 
the  following  differential  algebraic  form: 


dX  s 

Tt=f(x,y,u) 
0  =  h(x,y,  u ) 


(32) 
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Fig.  6.  Current/temperature  controller  configuration. 


where z  =  \j  a ^  a^]T  and  the  function  h(x, z,  u)  contains  Eqs.  (11), 
(19)  and  (20). 

4.  Closed-loop  dynamics  of  membrane  hydration 

Given  the  above  model,  we  can  now  proceed  to  analyze  the 
dynamic  behavior  of  the  membrane.  This  section  will  start  with 
a  very  simple  control-loop  structure  and  progressively  add  com¬ 
plexity.  The  purpose  of  this  progression  is  to  illustrate  the  coupling 
between  of  the  various  phenomena  within  the  fuel  cell  and  how 
these  are  impacted  by  the  various  levels  within  the  final  control 
loop  structure. 

4.1.  Current  control 

Under  this  scenario  the  only  form  of  regulation  will  be  with 
respect  to  current,  which  is  manipulated  by  changes  in  the  cell  (or 
load)  voltage.  The  configuration  of  Fig.  4  is  typical  of  the  electronic 


Fig.  8.  Expanded  view  of  Fig.  7. 


load  frequently  used  in  experimental  studies.  As  such,  we  have 
tuned  the  Proportional-Integral  (PI)  controller  to  be  fast  respond¬ 
ing. 

The  plots  of  Fig.  5  illustrate  the  response  to  step  changes  in  the 
current  density  set-point,  j(sp\  As  indicated  by  the  electrochemistry, 
a  decrease  in  cell  voltage  is  required  to  realize  the  desired  increase 
in  current.  The  increase  in  current  density  (and  thus  power  out¬ 
put)  will  increase  the  heat  production  rate,  Qgen,  as  observed  in  the 
temperature  plot  of  Fig.  5.  The  rise  in  solid  temperature  will  dramat¬ 
ically  impact  the  water  concentration  at  the  membrane  interfaces, 
^h2o  s  anc*  ch2o  s’  due  1:0  stron§  dependence  on  vapor  pressure. 


Fig.  7.  Simulation  of  current/temperature  controller. 
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Fig.  9.  Feed-forward  controller  configuration. 


This  change  in  surface  concentration  will  increase  the  flux  of  water 
from  the  membrane  to  the  cathode  gas  chamber.  This  is  in  con¬ 
trast  to  the  anode  side,  where  the  low  flow  rate  through  the  anode 
chamber  causes  the  moisture  content  of  the  anode  gas  to  track  the 
surface  concentration  and  results  in  a  nearly  uniform  flux  from  the 
anode  gas.  The  net  result  is  an  eventual  drying  out  of  the  membrane. 
The  impact  of  increased  water  production  and  electro-osmotic  drag 
can  be  seen  just  after  each  step  change.  This  is  observed  as  a  quick 
rise  in  water  content  at  the  cathode  interface,  A(rm),  as  well  as  a 
smaller  drop  at  the  anode,  A.(0).  The  eventual  drying  of  the  mem¬ 
brane  causes  its  resistance  to  increase,  which  the  current  controller 
compensates  for  by  dropping  cell  voltage.  However,  toward  the  end 
of  the  of  the  simulation  the  decrease  is  such  that  desired  current 
density  cannot  be  maintained. 


time  (sec) 


time  (sec) 


This  first  simulation  clearly  indicates  a  need  for  tempera¬ 
ture  control.  Although  temperature  control  is  expected  to  be 
part  of  any  fuel  cell  installation,  such  a  controller  is  likely  only 
able  to  regulate  the  bulk  (or  average)  temperature  of  the  stack. 
Given  the  spacial  nature  of  an  actual  fuel  cell  stack,  one  would 
expect  the  existence  of  local  hot-spots.  While  the  current  model 
cannot  capture  these  spacial  aspects,  the  above  simulation  sug¬ 
gests  the  type  of  phenomena  that  are  likely  occurring  at  the 
hot-spots. 

4.2.  Temperature  control 

We  now  consider  the  configuration  of  Fig.  6.  Under  this  sce¬ 
nario,  the  temperature  of  the  solid  is  regulated  by  manipulation 
of  the  cooling  jacket  flow.  Similar  to  the  previous  simulation,  a 
change  in  set-point  current  density  is  tracked  by  decreasing  cell 
voltage.  The  increase  in  power  output,  again,  result  in  an  increase 
in  heat  production.  However,  the  temperature  controller  responds 
by  increasing  jacket  flow  which  brings  the  cell  temperature  back 
to  the  set-point  (80  °C)  in  about  200  s.  Thus,  the  flux  of  water  from 
the  cathode  undergoes  a  much  smaller  increase.  This  coupled  with 
a  larger  increase  in  water  production  (due  to  a  larger  change  in  set- 
point  current  density)  results  in  an  eventual  increase  in  membrane 
hydration  (Fig.  7).  This  will  drop  ionic  resistance  and  allow  the  cur¬ 
rent  controller  to  achieve  its  set-point  at  a  higher  voltage.  Again  the 
electro-osmotic  drag  is  observed  at  the  step  changes.  At  the  anode 
interface,  A(0),  there  is  a  sudden  drop  in  hydration,  which  even- 


time  (sec) 


Fig.  10.  Simulation  with  feedforward  controller. 
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Fig.  11.  Simulation  with  feedforward  controller  and  Dm  =  5  x  10-6  cm  2  s-1 . 


tually  comes  back  due  to  the  overall  increase  in  hydration  level. 
At  the  cathode  interface,  A(rm),  we  see  more  structure  (see  Fig.  8). 
First  there  is  a  sudden  increase,  due  to  the  change  in  current  den¬ 
sity  and  thus  electro-osmotic  drag.  Then  the  profile  appears  to  level 
off,  due  to  an  increase  in  water  flux  to  the  cathode  gas.  Flowever,  as 
the  temperature  controller  kicks  in,  the  flux  is  brought  back  down 
and  the  hydration  level  again  increases,  which  again  increases  flux 
until  the  two  reach  a  new  equilibrium. 

4.3.  Flow  control 

In  the  previous  two  simulations,  the  flow  of  reactant  gases  to  the 
anode  and  cathode  gas  chambers  remained  unchanged  throughout. 
This,  however,  is  an  atypical  mode  of  operation.  The  more  com¬ 
mon  approach  is  to  vary  these  flowrates  based  on  a  fixed  reaction 
stoichiometry.  This  scheme  is  illustrated  by  the  feedforward  config¬ 
uration  of  Fig.  9.  Specifically,  the  inlet  flows  are  set  such  thatf£  =  2  * 
jWAmem/nFC“2  0  =  5pPkm3  s~3  and  F0C  =  2  */s^Amem/n^q :h  0  = 

40j(sP)cm3  s-3.  Fig.  10  illustrates  operation  under  this  scenario.  The 
impact  of  the  feedforward  action  is  observed  at  the  step  change 
times,  where  a  sharp  drop  in  anode  gas  temperature  and  water  con¬ 
centration  occurs.  Within  the  membrane  the  water  concentration 
at  the  anode  interface,  A(0),  again  begins  to  rise,  due  to  electro- 
osmotic  drag  and  increased  water  production.  But  then,  similar  to 
the  previous  case,  this  rise  is  cut  short  by  the  increase  in  water  flux 
from  the  membrane  due  to  the  rise  in  solid  temperature.  Flowever, 


in  contrast  to  the  previous  case,  the  flux  rise  is  greater  due  to  the 
drop  in  water  concentration  in  the  anode  gas  and  results  in  a  slight 
dip  in  membrane  hydration  at  the  anode  interface,  just  before  the 
eventual  rise  to  the  new  steady  state.  While  this  inverse  response 
is  quite  small  at  the  time  of  the  first  step  change  it  is  much  more 
pronounced  at  the  second  and  especially  the  third.  It  is  also  noted 
that  the  drop  in  anode  gas  water  concentration  results  in  lower 
hydration  levels  at  steady-state. 

5.  Discussion 

5.1.  Membrane  flooding 

The  last  two  simulations  (Figs.  7  and  10)  also  illustrate  an 
approach  to  the  flooding  condition.  Specifically,  the  membrane 
hydration  level  at  the  cathode,  X(rmem ),  approaches  the  critical  level 
of  A  =  14.  According  to  Eq.  (20),  such  a  value  of  A  will  result  in  a  water 
activity  in  the  membrane  of  Fig.  1  and  cause  the  partial  pressure  of 
water  at  the  surface  to  be  equal  the  vapor  pressure.  The  net  result 
is  a  saturation  in  the  water  removal  rate  in  the  form  of  vapor.  This 
suggests  that  the  only  mechanism  to  increase  the  rate  of  water 
removal  is  in  the  form  of  liquid.  It  should  additionally  be  noted 
that  A  <  14  does  not  preclude  the  existence  of  a  liquid  water  flux. 
Unfortunately,  the  literature  suggests  that  there  is  no  agreed  upon 
mechanism  describing  the  flux  of  liquid  water  from  the  membrane. 
As  such  the  current  model  makes  no  attempt  to  capture  the  flood- 
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ing  phenomenon.  We  should  also  note  that  temperature  gradients 
in  the  GDL  (specifically  a  cooling  near  the  GDL/gas  chamber  inter¬ 
face  )  could  cause  water  vapor  to  condense.  Again,  the  current  model 
makes  no  attempt  to  capture  this  phenomenon.  However,  this  con¬ 
densation  scenario  does  point  to  a  chain  of  events  that  would  lead 
to  flooding  within  the  spacially  dependent  realm  of  an  actual  fuel 
cell.  If  for  some  reason  a  cold  spot  were  to  occur  and  cause  vapor 
condensation  locally,  then  the  mass  transfer  of  the  oxygen  reactant 
would  also  be  reduced  locally.  This  reduction  in  reaction  rate  would 
cause  a  local  drop  in  heat  production  and  thus  further  local  cool¬ 
ing.  This  positive  feedback  would  continue  until  the  liquid  water 
grows  to  a  droplet  (and  eventually  a  slug)  and  is  moved  away  by 
the  hydrodynamics  of  the  flowing  gases.  While  the  current  model 
cannot  capture  this  chain  of  events,  the  above  discussion  does  point 
to  the  set  of  conditions  that  would  lead  to  flooding,  and  suggests 
the  future  development  of  a  predictive  type  controller  intended  to 
avoid  these  conditions  during  operation. 

5.2.  Impact  of  diffusion  in  membrane 

The  diffusion  coefficient  of  water  within  the  membrane,  Dm, 
plays  an  important  role  in  the  above  model.  In  the  previous  sim¬ 
ulations  we  used  the  Dm  value  suggested  by  O’Hayre  et  al.  [1] 
1.5  x  10-6  cm2  s-1.  However,  the  literature  [32]  suggests  a  lack  of 
agreement  on  the  value  of  this  parameter.  To  illustrate  sensitivity 
with  respect  to  Dm,  the  simulation  of  Fig.  1 0  was  repeated  (Fig.  1 1 ) 
using  a  different  value  for  Dm,  5  x  10-6  cm2  s-1.  The  curves  of  the 
second  simulation  have  nearly  identical  structure  to  those  of  the 
original  simulation.  Furthermore,  with  the  exception  of  membrane 
and  anode  chamber  water  concentration,  the  plot  values  are  about 
the  same.  The  main  difference  between  the  two  is  the  slope  of  the 
hydration  profile  in  the  membrane,  where  a  greater  diffusion  flux 
serves  to  flatten  the  profile.  With  regard  to  cathode  side  flooding, 
this  is  a  positive  outcome.  On  the  anode  side,  the  resulting  increase 
in  membrane  hydration  serves  to  reduce  the  flux  of  water  from 
the  anode  gas  chamber,  which  is  observed  as  an  increase  in  water 
in  the  anode  gas  chamber.  While  a  similar  consequence  should  be 
observed  in  the  cathode  chamber  (especially  due  to  the  drop  in 
membrane  hydration  at  the  cathode  interface),  the  impact  is  much 
less  pronounced,  due  to  the  larger  volumetric  flowrate  through  the 
cathode  gas  chamber. 

6.  Conclusions 

In  this  work  a  PEMFC  model  featuring  an  accumulation  based 
through  the  plane  membrane  characterization  was  combined  with 
a  variety  closed-loop  control  structures.  The  simulations  presented 
illustrate  a  complex  set  of  possible  responses,  owing  to  the  interac¬ 
tion  of  multiple  phenomena  (electro-chemical,  chemical,  thermal, 
and  membrane  hydration)  occurring  at  multiple  time-scales.  The 
addition  of  regulatory  and  feedforward  control  loops  was  shown  to 


have  a  significant  impact  on  response  structure  and  settling-time, 
and  as  such  should  be  included  in  the  dynamic  characterization  of 
the  PEMFC.  It  was  also  noted  that  changes  in  membrane  diffusivity 
had  a  significant  impact  on  water  accumulation  levels  within  the 
membrane. 
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